In this Notebook I want to explore methods of phylogenetic imputation, that is, using phylogenetic information to impute missing values in any given dataset, be these continuous (including count) or discrete.
It is important to understand that most methods implicitly rely on the assumption that the phylogeny will help predict the traits better. Either because closely related species tend to share more similar trait values than distantly related species, or because we assume that traits are constrained by evolution according to a certain model (like Brownian Motion). This is of course not always true. Thus is it also important to test for phylogenetic signal in any given trait (more on this later). If traits do not show any phylogenetic signal, then trying to impute them using the phylogeny becomes uninformative and can potentially reduce the precision of the imputed values.
One personal observation is that phylogenetic signal will increase with the size of the dataset. That is, if you are trying to predict traits for a list of 100 species, you might find that using a phylogeny with only 100 tips will perform poorly relative to a phylogeny of 1,000 tips.
It is also possible to incorporate a correlation structure linked to other traits that might be helpful in imputing the missing values, just like you would do if you used a linear model, or imputation by PCA or even Multiple Imputation by Chained Equations (MICE).
Finally, it will be important to have a measure of validation of the imputed values. A cross validation or monte carlo procedure can be used for this (see Molina-Venegas, 2023).
Load the libraries
library(phytools)
library(phangorn)
library(Rphylopars)
library(tidyverse)
library(ape)
library(geiger)
library(missForest)
library(funspace)
library(readxl)
library(car)
library(picante)
library(Hmisc)
library(broom)
library(PerformanceAnalytics)
library(ade4)
Ok, let’s get started with some simple examples.
Lets assume that our dataset consists of continuous and discrete data. Here we can start by using the dataset from Sandra diaz et al., 2022 which is an enhanced dataset containing the 6 major traits of the global spectrum of plant form and function (Diaz et al. 2016) plus some other information for about 46,000 species.
data <- read_xlsx("../Functional traits/Diaz et al. 2022. Enhanced Trait database/Try2023123184331480_Dataset/Dataset/Species_mean_traits.xlsx", sheet = 1, col_names = T) %>%
dplyr::rename(Species = `Species name standardized against TPL`,
LDMC = `LDMC (g/g)`,
LA=`Leaf area (mm2)`,
Nmass=`Nmass (mg/g)`,
LMA = `LMA (g/m2)`,
Height = `Plant height (m)`,
Seedmass = `Diaspore mass (mg)`,
SSD = `SSD combined (mg/mm3)`)
data$Species <- gsub(" ", "_", data$Species)
data <- strings2factors(data)
The following character variables were converted to factors
Taxonomic level Status according to TPL Genus Family Phylogenetic Group within angiosperms Phylogenetic Group General Adaptation to terrestrial or aquatic habitats Woodiness Growth Form Nutrition type (parasitism) Leaf type
str(data)
tibble [46,047 Ă— 32] (S3: tbl_df/tbl/data.frame)
$ TRY 30 AccSpecies ID : num [1:46047] 89902 2 42370 4 73137 ...
$ Species : chr [1:46047] "Aaronsohnia_pubescens" "Abarema_adenophora" "Abarema_adenophorum" "Abarema_brachystachya" ...
$ Taxonomic level : Factor w/ 3 levels "genus","species",..: 2 2 2 2 2 2 2 2 2 2 ...
$ Status according to TPL : Factor w/ 5 levels "accepted","invalid name",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Genus : Factor w/ 7227 levels "Aaronsohnia",..: 1 2 2 2 2 2 2 2 2 2 ...
$ Family : Factor w/ 419 levels "Acanthaceae",..: 100 214 214 214 214 214 214 214 214 214 ...
$ Phylogenetic Group within angiosperms : Factor w/ 4 levels "Angiosperm_Magnoliid",..: 3 3 3 3 3 3 3 3 3 3 ...
$ Phylogenetic Group General : Factor w/ 3 levels "Angiosperm","Gymnosperm",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Adaptation to terrestrial or aquatic habitats: Factor w/ 4 levels "aquatic","aquatic/semiaquatic",..: 4 4 4 4 4 4 4 4 4 4 ...
$ Woodiness : Factor w/ 4 levels "non-woody","semi-woody",..: 1 3 3 3 3 3 3 3 3 3 ...
$ Growth Form : Factor w/ 11 levels "bamboo graminoid",..: 5 11 NA 9 11 11 11 11 11 11 ...
$ Succulence : logi [1:46047] NA NA NA NA NA NA ...
$ Nutrition type (parasitism) : Factor w/ 4 levels "hemiparasitic",..: 3 3 NA 3 3 3 3 3 3 3 ...
$ Nutrition type (carnivory) : logi [1:46047] NA NA NA NA NA NA ...
$ Leaf type : Factor w/ 5 levels "broadleaved",..: 1 1 1 1 NA NA 1 1 1 1 ...
$ LA : num [1:46047] NA 7162 NA NA 700 ...
$ Leaf area (n.o.) : num [1:46047] NA 5 NA NA 1 1 16 NA NA NA ...
$ Nmass : num [1:46047] NA 24.4 NA 30.4 NA ...
$ Nmass (n.o.) : num [1:46047] NA 6 NA 1 NA NA 18 2 1 22 ...
$ LMA : num [1:46047] NA 87.1 NA NA NA ...
$ LMA (n.o.) : num [1:46047] NA 4 NA NA NA NA 14 2 6 7 ...
$ Height : num [1:46047] 0.2 NA NA 10 NA ...
$ Plant height (n.o.) : num [1:46047] 2 NA NA 1 NA NA 13 NA 2 7 ...
$ Seedmass : num [1:46047] NA NA NA NA NA ...
$ Diaspore mass (n.o.) : num [1:46047] NA NA NA NA NA NA 16 NA 7 NA ...
$ SSD observed (mg/mm3) : num [1:46047] NA 0.648 0.36 0.53 NA ...
$ SSD (n.o.) : num [1:46047] NA 4 2 2 NA NA 29 NA 3 7 ...
$ LDMC : num [1:46047] NA NA NA NA NA NA NA NA NA NA ...
$ LDMC (n.o.) : num [1:46047] NA NA NA NA NA NA NA NA NA NA ...
$ SSD imputed (mg/mm3) : num [1:46047] NA NA NA NA NA NA NA NA NA NA ...
$ SSD : num [1:46047] NA 0.648 0.36 0.53 NA ...
$ Number of traits with values : num [1:46047] 1 4 1 3 1 1 6 2 5 4 ...
We can start by filtering this datset to select only the species which have all observations of traits and are subgeneric level observations.
filtered_data <- data %>% filter(`Number of traits with values` > 5 & `Taxonomic level` != "genus")
dim(filtered_data)
[1] 2214 32
Ok so that’s 2214 species left.
Check if there are any NAs
cat(sum(is.na(filtered_data)) / (nrow(filtered_data)*ncol(filtered_data)) * 100, "%")
15.74498 %
So still about 15% missing data, that is from the column LDMC.
Check the distribution of each variable.
filtered_data %>% dplyr::select(LDMC, LA, Nmass, LMA, SSD, Height, Seedmass) %>%
gather %>%
ggplot(aes(value, fill=key)) + geom_density() + facet_wrap(facets = ~key, scales = "free")
Now log them all and check they are more or less normally distributed.
new <- filtered_data %>% column_to_rownames("Species") %>%
dplyr::select(LDMC, LA, Nmass, LMA, SSD, Height, Seedmass) %>% log()
new %>%
gather %>%
ggplot(aes(value, fill=key)) + geom_density() + facet_wrap(facets = ~key, scales = "free")
# Check for normality
apply(new, 2, shapiro.test) %>% do.call(rbind, .)
statistic p.value method data.name
LDMC 0.9933977 0.0004499156 "Shapiro-Wilk normality test" "newX[, i]"
LA 0.9829842 1.313264e-15 "Shapiro-Wilk normality test" "newX[, i]"
Nmass 0.9924938 2.978751e-09 "Shapiro-Wilk normality test" "newX[, i]"
LMA 0.9864945 1.262429e-13 "Shapiro-Wilk normality test" "newX[, i]"
SSD 0.9598879 2.850217e-24 "Shapiro-Wilk normality test" "newX[, i]"
Height 0.9232025 3.320459e-32 "Shapiro-Wilk normality test" "newX[, i]"
Seedmass 0.9735828 8.540292e-20 "Shapiro-Wilk normality test" "newX[, i]"
Ok so it’s not perfect but it is much better.
Let’s see what the correlations between traits are.
new %>%
as.matrix() %>%
Hmisc::rcorr() %>%
broom::tidy()
new %>%
chart.Correlation()
So there is some good correlation between log(SSD) and log(LDMC) and between log(SSD) and log(Height). log(Height) and log(Seedmass) aren’ too bad either.
Ok now we need a phylogeny for these. Fortunately for us, the package funspace includes a phylogeny for these species Carmona et al. 2021.
phy <- funspace::phylo
Most tests for phylogenetic signal rely on Brownian Motion (BM),
which is equivalent to a random walk through time with a known variance
sigma^2. This means that as time increases, the variance
also increases between two given tips.
Two widely used metrics for measuring phylogenetic signal under BM
are Pagel’s \(\lambda\) and Bloomberg’s
K (see Pagel 1999 and Bloomberg 2003). These are usually
defined between [0,1] such as 0 is complete randomness and 1 perfectly
matches BM. \(\lambda\) can effectively
be higher than 1 but is generally no well defined over that limit.
K can be higher than 1, effectively meaning there is
stronger signal than you would expect under BM.
For binary categorical traits, the only method I am confident in
using is the phylo.d method by Fritz and Purvis (2010), implemented in
caper. Another method is the \(\delta\) statistic by Borges et al. 2019
(implemented here https://github.com/mrborges23/delta_statistic/tree/master).
Let’s see if any of our traits can be measured with these two metrics.
As a reminder, I am using log transformed variables for modelling. This is easier to fit under Brownian models of evolution and so is more appropriate to test for phylogenetic signal under these models.
First, let’s match our phylogeny and trait data. Because the dataset is quite large, let’s subsample 500 species and run our analyses on this. We can expand the dataset to see how the tests perform on a larger dataset.
Quick fix in the tree tip labels
phy$tip.label[4884] <- "Cercocarpus_montanus_var._glaber"
Ok now we can compute K and \(\lambda\) for each trait. Some models work
best when the trees are completely resolved (ie. no polytomies). Thus we
will make sure our tree fits this assumption. If necessary we can
transform the tree with polytomies into dichotomies with zero branch
lengths using the hand multi2di function.
is.binary(phylo_traits$phy)
[1] FALSE
phylo_traits$phy <- multi2di(phylo_traits$phy)
Great, now let’s calculate K and \(\lambda\) for Leaf Mass per Area.
LMA <- phylo_traits$data$LMA %>% setNames(rownames(phylo_traits$data))
(K.LMA<-phylosig(tree = phylo_traits$phy, x = LMA, method = "K", test = T, nsim=500) )
Phylogenetic signal K : 0.0761484
P-value (based on 500 randomizations) : 0.002
plot.phylosig(K.LMA)
# This takes a very long time to compute for larger datasets but the geiger function `fitcontinuous` with model = lambda, provides the same result
(lambda.LMA<-phylosig(tree = phylo_traits$phy, x = LMA, method = "lambda", test = T) )
Phylogenetic signal lambda : 0.859225
logL(lambda) : -378.058
LR(lambda=0) : 124.506
P-value (based on LR test) : 6.5287e-29
Ok, we can see that there is a very low value of K, and that this value is significantly different from random.
But we can also see that there is a much higher value for \(\lambda\).
How can this be? Both are supposed to quantify phylogenetic signal. They do in fact look at different aspects of phylogenetic signal. K is a ratio of within to between clade variance, whilst lambda is a transformation of the diagonal element of the pVCV matrix.
The p-values here compare the results to a null expectation of randomness, but there is not automated way of comparing this to a null expectation of Brownian motion. We have to do this manually.
#p.93 of the book Phylogenetic Comparative Methods in R (Revell & Harmon, 2022)
#simulate 1000 datasets
nullX <- fastBM(phylo_traits$phy, nsim = 1000)
## for each, carry out a test of phylogenetic signal
## amd accumulate these into a vector using apply
nullK <- apply(nullX,2,phylosig, tree=phylo_traits$phy)
#calcultae p.values
Pvals_LMA <- mean(nullK<=K.LMA$K)
Pvals_LMA
[1] 0
Check this out visually
hist(c(nullK, K.LMA$K), breaks=100, col="lightgray", border="black", main="", xlab="K", las=1, cex.axis=0.7, cex.lab=0.9, ylim=c(0,2000), xlim=c(0,4))
arrows(x0=K.LMA$K, y0=par()$usr[4], y1=0, length=0.12, col=make.transparent("blue", 0.5), lwd=2)
text(K.LMA$K, 0.96*par()$usr[4], paste("observed value of K (P= ", round(Pvals_LMA, 4), ")", sep = ""), pos=4, cex=0.8)
We can see that the distribution of K under Brownian motion is very widespread! But our value of K is indeed lower than expected under Brownian motion.
What about the value of lambda? Can we reject a null of lambda = 1?
We can use a likelihood ratio test for this one.
(LR_LMA <- -2*(lambda.LMA$lik(1) - lambda.LMA$logL) )
[1] 469.5403
Now compute the p-value (assuming chi-square distribution under the null hypothesis of lambda = 1)
Pval_lambda_LMA <- pchisq(LR_LMA, df=1, lower.tail = F )
Pval_lambda_LMA
[1] 4.034683e-104
So we can reject the null that \(\lambda\) = 1.
This means that both K and \(\lambda\) are telling us there is more
phylogenetic signal than random, but less than expected under Brownian
motion. So perhaps overdispersed ?
We can also compare these results to other models of evolution, like Early-Burst (EB) or Ornstein-Uhlenbeck (OU) models.
(BM.LMA <- fitContinuous(phylo_traits$phy, LMA, model = "BM") )
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.020824
z0 = 4.520216
model summary:
log-likelihood = -612.828125
AIC = 1229.656250
AICc = 1229.680395
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.LMA <- fitContinuous(phylo_traits$phy, LMA, model = "EB") )
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.020832
z0 = 4.520218
model summary:
log-likelihood = -612.832999
AIC = 1231.665999
AICc = 1231.714386
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 1
frequency of best fit = 0.010
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.LMA <- fitContinuous(phylo_traits$phy, LMA, model = "OU") )
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 0.145550
sigsq = 0.105084
z0 = 4.199499
model summary:
log-likelihood = -433.949344
AIC = 873.898688
AICc = 873.947075
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 44
frequency of best fit = 0.440
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(Lambda.LMA <- fitContinuous(phylo_traits$phy, LMA, model = "lambda") )
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.859240
sigsq = 0.002478
z0 = 4.514549
model summary:
log-likelihood = -378.057990
AIC = 762.115979
AICc = 762.164366
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 16
frequency of best fit = 0.160
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
It seems like the a parameter for the EB model is close
to 0, which is effectively equivalent to a BM.
If we check the output for the lambda model, this is the same as
using the phylosig function from phytools but does it much
quicker and with a few more bits of information.
#However if the OU model seems to be at bounds with the `alpha` parameter, then we can increase this and see how the model performs.
#(OU.LMA <- fitContinuous(phylo_traits$phy, LMA, model = "OU", bounds = list(alpha=c(0,50))) )
It’s also good to see how these compare to random noise model.
white is a white-noise (non-phylogenetic) model, which
assumes data come from a single normal distribution with no covariance
structure among species. The variance parameter sigsq takes the same
bounds defined under the BM model.
null.LMA <- fitContinuous(phy = phylo_traits$phy, dat = LMA, model = "white")
Now we can compare how well each model performs using AIC
(aics_models <- setNames(c(AIC(BM.LMA), AIC(EB.LMA), AIC(OU.LMA), AIC(null.LMA)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
1229.6562 1231.6660 873.8987 884.6218
Seems like the OU model has the lowest AIC. We can also check the weights.
aic.w(aics_models)
BM EB OU NULL
0.00000000 0.00000000 0.99532836 0.00467164
What about comparing OU to Lambda model ?
aic.w(setNames(c(AIC(OU.LMA), AIC(Lambda.LMA)), c("OU", "Lambda")) )
OU Lambda
0 1
So in this case the lambda transformation fits the data better.
We can check this out visually and apply the lambda and OU transformations.
#normal tree
contMap(tree = phylo_traits$phy, x = LMA, ftype = "off", type = "fan", lwd=2)
#lambda transformation
#contMap(tree = lambdaTree(phylo_traits$phy,lambda = Lambda.LMA$opt$lambda), x = LMA, fsize = 0.7, ftype = "off", type = "fan", lwd=2)
# or use this
contMap(tree = rescale(phylo_traits$phy, model = "lambda", lambda = Lambda.LMA$opt$lambda), x = LMA, ftype = "off", type = "fan", lwd=2)
#OU transformation
contMap(tree = rescale(phylo_traits$phy, model = "OU", alpha=OU.LMA$opt$alpha), x = LMA, ftype = "off", type = "fan", lwd=2)
The OU model seems quite strange. It looks like most branches have been shortened so much they are almost like a star phylogeny.
We can do this for all the traits one by one and see how much phylogenetic signal they have.
SSD <- phylo_traits$data$SSD %>% setNames(rownames(phylo_traits$data))
Height <- phylo_traits$data$Height %>% setNames(rownames(phylo_traits$data))
LA <- phylo_traits$data$LA%>% setNames(rownames(phylo_traits$data))
LDMC <- phylo_traits$data$LDMC %>% setNames(rownames(phylo_traits$data)) # has some NA values
Nmass <- phylo_traits$data$Nmass %>% setNames(rownames(phylo_traits$data))
Seedmass <- phylo_traits$data$Seedmass %>% setNames(rownames(phylo_traits$data))
Now fit the models
(K.SSD<-phylosig(tree = phylo_traits$phy, x = SSD, method = "K", test = T, nsim=500) )
Phylogenetic signal K : 0.134274
P-value (based on 500 randomizations) : 0.002
plot.phylosig(K.SSD)
(Lambda.SSD <- fitContinuous(phylo_traits$phy, SSD, model = "lambda") )
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.957333
sigsq = 0.002486
z0 = -1.245813
model summary:
log-likelihood = -242.942713
AIC = 491.885426
AICc = 491.933814
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 1
frequency of best fit = 0.010
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(BM.SSD <- fitContinuous(phylo_traits$phy, SSD, model = "BM") )
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.008738
z0 = -1.250112
model summary:
log-likelihood = -395.732158
AIC = 795.464317
AICc = 795.488462
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.SSD <- fitContinuous(phylo_traits$phy, SSD, model = "EB") )
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.008742
z0 = -1.250115
model summary:
log-likelihood = -395.736140
AIC = 797.472280
AICc = 797.520667
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 1
frequency of best fit = 0.010
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.SSD <- fitContinuous(phylo_traits$phy, SSD, model = "OU") )
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 0.019912
sigsq = 0.014408
z0 = -1.096850
model summary:
log-likelihood = -328.794977
AIC = 663.589953
AICc = 663.638340
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 19
frequency of best fit = 0.190
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(null.SSD <- fitContinuous(phylo_traits$phy, SSD, model = "white"))
GEIGER-fitted comparative model of continuous data
fitted ‘white’ model parameters:
sigsq = 0.320221
z0 = -1.124992
model summary:
log-likelihood = -424.783457
AIC = 853.566915
AICc = 853.591060
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
Here \(\lambda\) is 0.96, so pretty strong
Check which fits best and the AIC weights
(aics_models <- setNames(c(AIC(BM.SSD), AIC(EB.SSD), AIC(OU.SSD), AIC(null.SSD)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
795.4643 797.4723 663.5900 853.5669
aic.w(aics_models)
BM EB OU NULL
0 0 1 0
Again, OU seems to do best.
(K.Nmass<-phylosig(tree = phylo_traits$phy, x = Nmass, method = "K", test = T, nsim=500) )
Phylogenetic signal K : 0.0108426
P-value (based on 500 randomizations) : 0.554
plot.phylosig(K.Nmass)
(Lambda.Nmass <- fitContinuous(phylo_traits$phy, Nmass, model = "lambda") )
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.800613
sigsq = 0.000846
z0 = 2.849721
model summary:
log-likelihood = -159.493136
AIC = 324.986273
AICc = 325.034660
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 14
frequency of best fit = 0.140
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(BM.Nmass <- fitContinuous(phylo_traits$phy, Nmass, model = "BM") )
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.058183
z0 = 2.853405
model summary:
log-likelihood = -869.700663
AIC = 1743.401326
AICc = 1743.425471
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.Nmass <- fitContinuous(phylo_traits$phy, Nmass, model = "EB") )
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.058206
z0 = 2.853405
model summary:
log-likelihood = -869.706608
AIC = 1745.413216
AICc = 1745.461603
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 1
frequency of best fit = 0.010
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.Nmass <- fitContinuous(phylo_traits$phy, Nmass, model = "OU") )
Warning:
Parameter estimates appear at bounds:
alpha
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 2.718282
sigsq = 0.861545
z0 = 3.027972
model summary:
log-likelihood = -247.247537
AIC = 500.495075
AICc = 500.543462
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 38
frequency of best fit = 0.380
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(null.Nmass <- fitContinuous(phylo_traits$phy, Nmass, model = "white"))
GEIGER-fitted comparative model of continuous data
fitted ‘white’ model parameters:
sigsq = 0.149826
z0 = 3.028794
model summary:
log-likelihood = -234.898773
AIC = 473.797546
AICc = 473.821691
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
Here K seems to be non-significant. Lambda is moderate to high.
Check AIC weights
(aics_models <- setNames(c(AIC(BM.Nmass), AIC(EB.Nmass), AIC(OU.Nmass), AIC(null.Nmass)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
1743.4013 1745.4132 500.4951 473.7975
aic.w(aics_models)
BM EB OU NULL
0.00000000 0.00000000 0.00000159 0.99999841
This time a null model of white noise (grand mean) is best.
(K.Height<-phylosig(tree = phylo_traits$phy, x = Height, method = "K", test = T, nsim=500) )
Phylogenetic signal K : 0.198507
P-value (based on 500 randomizations) : 0.002
plot.phylosig(K.Height)
(Lambda.Height <- fitContinuous(phylo_traits$phy, Height, model = "lambda") )
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.982437
sigsq = 0.033418
z0 = 0.807796
model summary:
log-likelihood = -826.258043
AIC = 1658.516086
AICc = 1658.564473
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 12
frequency of best fit = 0.120
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(BM.Height <- fitContinuous(phylo_traits$phy, Height, model = "BM") )
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.079894
z0 = 0.807071
model summary:
log-likelihood = -948.979485
AIC = 1901.958971
AICc = 1901.983116
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.Height <- fitContinuous(phylo_traits$phy, Height, model = "EB") )
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.079927
z0 = 0.807065
model summary:
log-likelihood = -948.982876
AIC = 1903.965752
AICc = 1904.014139
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 1
frequency of best fit = 0.010
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.Height <- fitContinuous(phylo_traits$phy, Height, model = "OU") )
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 0.012628
sigsq = 0.115850
z0 = 0.837843
model summary:
log-likelihood = -908.809345
AIC = 1823.618690
AICc = 1823.667077
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 9
frequency of best fit = 0.090
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(null.Height <- fitContinuous(phylo_traits$phy, Height, model = "white"))
GEIGER-fitted comparative model of continuous data
fitted ‘white’ model parameters:
sigsq = 4.373830
z0 = 0.399501
model summary:
log-likelihood = -1078.379066
AIC = 2160.758132
AICc = 2160.782277
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
Lambda is pretty much equal to 1
Check AIC weights
(aics_models <- setNames(c(AIC(BM.Height), AIC(EB.Height), AIC(OU.Height), AIC(null.Height)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
1901.959 1903.966 1823.619 2160.758
aic.w(aics_models)
BM EB OU NULL
0 0 1 0
OU model fits data best
(K.LA<-phylosig(tree = phylo_traits$phy, x = LA, method = "K", test = T, nsim=500) )
Phylogenetic signal K : 0.128972
P-value (based on 500 randomizations) : 0.002
plot.phylosig(K.LA)
(Lambda.LA <- fitContinuous(phylo_traits$phy, LA, model = "lambda") )
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.875074
sigsq = 0.031594
z0 = 5.328225
model summary:
log-likelihood = -998.321917
AIC = 2002.643833
AICc = 2002.692220
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 22
frequency of best fit = 0.220
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(BM.LA <- fitContinuous(phylo_traits$phy, LA, model = "BM") )
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.171152
z0 = 5.340679
model summary:
log-likelihood = -1139.441718
AIC = 2282.883436
AICc = 2282.907581
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.LA <- fitContinuous(phylo_traits$phy, LA, model = "EB") )
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.171218
z0 = 5.340664
model summary:
log-likelihood = -1139.446085
AIC = 2284.892170
AICc = 2284.940557
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 13
frequency of best fit = 0.130
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.LA <- fitContinuous(phylo_traits$phy, LA, model = "OU") )
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 0.040244
sigsq = 0.388378
z0 = 6.811381
model summary:
log-likelihood = -1033.938724
AIC = 2073.877447
AICc = 2073.925834
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 20
frequency of best fit = 0.200
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(null.LA <- fitContinuous(phylo_traits$phy, LA, model = "white"))
GEIGER-fitted comparative model of continuous data
fitted ‘white’ model parameters:
sigsq = 4.365239
z0 = 6.738389
model summary:
log-likelihood = -1077.887512
AIC = 2159.775023
AICc = 2159.799168
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
Check AIC weights
(aics_models <- setNames(c(AIC(BM.LA), AIC(EB.LA), AIC(OU.LA), AIC(null.LA)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
2282.883 2284.892 2073.877 2159.775
aic.w(aics_models)
BM EB OU NULL
0 0 1 0
Again, OU
(K.LDMC<-phylosig(tree = phylo_traits$phy, x = LDMC[!is.na(LDMC)], method = "K", test = T, nsim=500) )
[1] "some species in tree are missing from x , dropping missing taxa from the tree"
Phylogenetic signal K : 0.216002
P-value (based on 500 randomizations) : 0.002
plot.phylosig(K.LDMC)
(Lambda.LDMC <- fitContinuous(phylo_traits$phy, LDMC[!is.na(LDMC)], model = "lambda") )
Warning: The following tips were not found in 'data' and were dropped from 'phy':
Abies_grandis
Acacia_aroma
Acacia_auriculiformis
Acacia_doratoxylon
Acacia_karroo
Acacia_mearnsii
Acacia_praecox
Acer_monspessulanum
Acer_saccharinum
Acer_saccharum
Albizia_saman
Alliaria_petiolata
Alnus_rhombifolia
Ambrosia_psilostachya
Anagallis_tenella
Andira_inermis
Anemone_nemorosa
Annona_spraguei
Apeiba_tibourbou
Araucaria_araucana
Arenaria_serpyllifolia
Artemisia_californica
Artemisia_dracunculus
Aspidosperma_album
Banksia_marginata
Betula_papyrifera
Betula_platyphylla
Bombax_ceiba
Briza_media
Bromus_hordeaceus
Bromus_sterilis
Brosimum_utile
Brownea_grandiceps
Butea_monosperma
Byrsonima_crassifolia
Calamagrostis_epigejos
Calophyllum_longifolium
Campanula_rotundifolia
Capsella_bursa-pastoris
Carex_binervis
Carex_hirta
Carex_pallescens
Carex_panicea
Carya_glabra
Caryocar_glabrum
Casearia_sylvestris
Castanea_dentata
Ceanothus_spinosus
Cecropia_obtusifolia
Cecropia_peltata
Ceiba_speciosa
Cercocarpus_montanus_var._glaber
Cinnamomum_subavenium
Cirsium_wheeleri
Cojoba_rufescens
Colophospermum_mopane
Cordia_bicolor
Cordia_caffra
Cordia_megalantha
Cornus_glabrata
Cornus_sericea
Crepis_pyrenaica
Croton_capitatus
Cupania_rufescens
Cupressus_sempervirens
Cussonia_spicata
Cytisus_scoparius
Dacrydium_cupressinum
Danthonia_decumbens
Daphniphyllum_pentandrum
Drimys_winteri
Dryandra_sessilis
Drypetes_standleyi
Duguetia_quitarensis
Ekebergia_capensis
Elaeagnus_rhamnoides
Emmotum_fagifolium
Enterolobium_schomburgkii
Epilobium_palustre
Equisetum_arvense
Equisetum_fluviatile
Equisetum_palustre
Erica_tetralix
Erythrophleum_chlorostachys
Eschweilera_decolorans
Eschweilera_parviflora
Eschweilera_sagotiana
Eucalyptus_acmenoides
Eucalyptus_grandis
Eucalyptus_umbra
Euphorbia_brachycera
Euphorbia_corollata
Faramea_occidentalis
Fragaria_vesca
Frangula_californica
Fraxinus_uhdei
Galium_verum
Garcinia_macrophylla
Grewia_monticola
Griselinia_littoralis
Guatteria_dumetorum
Gymnosporia_heterophylla
Hakea_cygna
Hakea_preissii
Hamamelis_virginiana
Hampea_appendiculata
Hampea_nutricia
Heliocarpus_appendiculatus
Heracleum_sphondylium
Heteromeles_arbutifolia
Holcus_mollis
Holodiscus_discolor
Ilex_mitis
Inga_acreana
Inga_goldmanii
Inga_nobilis
Inga_sapindoides
Iris_missouriensis
Iryanthera_juruensis
Jacaranda_copaia
Jacaratia_digitata
Jacobaea_vulgaris
Juncus_articulatus
Juniperus_communis
Kigelia_africana
Koeleria_macrantha
Lacmellea_panamensis
Lactuca_muralis
Laetia_corymbulosa
Larix_kaempferi
Larrea_cuneifolia
Larrea_divaricata
Lathyrus_pratensis
Lecythis_zabucajo
Leiospermum_racemosum
Lepechinia_calycina
Leptospermum_polygalifolium
Leucadendron_salignum
Ligustrum_lucidum
Liquidambar_styraciflua
Lotus_pedunculatus
Lysimachia_vulgaris
Lythrum_salicaria
Macaranga_hypoleuca
Machilus_chinensis
Macrolobium_gracile
Magnolia_grandiflora
Malus_sylvestris
Mangifera_indica
Margaritaria_nobilis
Melaleuca_leucadendra
Metrosideros_polymorpha
Mimulus_aurantiacus
Molinia_caerulea
Mosannona_garwoodii
Murraya_paniculata
Myrcia_splendens
Nothofagus_betuloides
Nothofagus_dombeyi
Nothofagus_fusca
Nothofagus_solandri
Ochroma_pyramidale
Ocotea_floribunda
Ocotea_puberula
Opuntia_sulphurea
Ormosia_macrocalyx
Orthion_oblanceolatum
Ouratea_lucens
Oxalis_acetosella
Pachira_aquatica
Pachira_insignis
Packera_multilobata
Palicourea_tetragona
Parinari_campestris
Parkinsonia_praecox
Parnassia_palustris
Peltogyne_venosa
Perebea_xanthochyma
Petrea_volubilis
Picea_sitchensis
Picramnia_latifolia
Pimpinella_major
Pinus_clausa
Pinus_koraiensis
Pinus_massoniana
Piper_aequale
Piscidia_carthagenensis
Pistacia_lentiscus
Pongamia_pinnata
Populus_nigra
Populus_tremuloides
Pourouma_bicolor
Pouteria_cladantha
Pouteria_ephedrantha
Pradosia_cochlearia
Prosopis_flexuosa
Prosopis_torquata
Protium_tenuifolium
Prunus_mahaleb
Prunus_pensylvanica
Pseudocymopterus_montanus
Psidium_guajava
Psychotria_marginata
Pterocarpus_rohrii
Pterygota_alata
Quararibea_asterolepis
Quassia_amara
Quercus_lusitanica
Quercus_robur
Quercus_shumardii
Quercus_stellata
Quercus_variabilis
Ranunculus_bulbosus
Reynoutria_japonica
Rhamnus_alaternus
Rollinia_pittieri
Rosa_woodsii
Rubus_parviflorus
Ruellia_humilis
Rumex_obtusifolius
Salix_cinerea
Salix_lasiolepis
Salix_nigra
Schinopsis_marginata
Schizolobium_parahyba
Sequoia_sempervirens
Shorea_leprosula
Shorea_macroptera
Shorea_parvifolia
Shorea_smithiana
Silphium_terebinthinaceum
Sloanea_guianensis
Sloanea_terniflora
Solanum_dulcamara
Sorbus_aria
Sorocea_steinbachii
Spondias_mombin
Spondias_pinnata
Sporobolus_interruptus
Sterculia_speciosa
Stipa_comata
Tabernaemontana_arborea
Tabernaemontana_donnell-smithii
Tachigali_versicolor
Talisia_princeps
Tanacetum_corymbosum
Tapirira_obtusa
Taraxacum_campylodes
Taxodium_distichum
Taxus_baccata
Terminalia_amazonia
Theobroma_cacao
Theobroma_speciosum
Thevetia_ahouai
Trattinnickia_aspera
Triadica_cochinchinensis
Trichilia_pleeana
Trichilia_quadrijuga
Trichospermum_mexicanum
Trifolium_dubium
Triplaris_cumingiana
Ulmus_glabra
Ulmus_minor
Ulmus_pumila
Ulmus_rubra
Unonopsis_rufescens
Vaccinium_myrtillus
Vaccinium_vitis-idaea
Verbena_bracteata
Veronica_chamaedrys
Viburnum_opulus
Victoria_amazonica
Viola_canina
Vismia_baccifera
Vitex_negundo
Vitex_pinnata
Ziziphus_mucronata
Zuelania_guidonia
GEIGER-fitted comparative model of continuous data
fitted ‘lambda’ model parameters:
lambda = 0.740265
sigsq = 0.001569
z0 = -1.707557
model summary:
log-likelihood = -80.434508
AIC = 166.869016
AICc = 166.979108
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 20
frequency of best fit = 0.200
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(BM.LDMC <- fitContinuous(phylo_traits$phy, LDMC[!is.na(LDMC)], model = "BM") )
Warning: The following tips were not found in 'data' and were dropped from 'phy':
Abies_grandis
Acacia_aroma
Acacia_auriculiformis
Acacia_doratoxylon
Acacia_karroo
Acacia_mearnsii
Acacia_praecox
Acer_monspessulanum
Acer_saccharinum
Acer_saccharum
Albizia_saman
Alliaria_petiolata
Alnus_rhombifolia
Ambrosia_psilostachya
Anagallis_tenella
Andira_inermis
Anemone_nemorosa
Annona_spraguei
Apeiba_tibourbou
Araucaria_araucana
Arenaria_serpyllifolia
Artemisia_californica
Artemisia_dracunculus
Aspidosperma_album
Banksia_marginata
Betula_papyrifera
Betula_platyphylla
Bombax_ceiba
Briza_media
Bromus_hordeaceus
Bromus_sterilis
Brosimum_utile
Brownea_grandiceps
Butea_monosperma
Byrsonima_crassifolia
Calamagrostis_epigejos
Calophyllum_longifolium
Campanula_rotundifolia
Capsella_bursa-pastoris
Carex_binervis
Carex_hirta
Carex_pallescens
Carex_panicea
Carya_glabra
Caryocar_glabrum
Casearia_sylvestris
Castanea_dentata
Ceanothus_spinosus
Cecropia_obtusifolia
Cecropia_peltata
Ceiba_speciosa
Cercocarpus_montanus_var._glaber
Cinnamomum_subavenium
Cirsium_wheeleri
Cojoba_rufescens
Colophospermum_mopane
Cordia_bicolor
Cordia_caffra
Cordia_megalantha
Cornus_glabrata
Cornus_sericea
Crepis_pyrenaica
Croton_capitatus
Cupania_rufescens
Cupressus_sempervirens
Cussonia_spicata
Cytisus_scoparius
Dacrydium_cupressinum
Danthonia_decumbens
Daphniphyllum_pentandrum
Drimys_winteri
Dryandra_sessilis
Drypetes_standleyi
Duguetia_quitarensis
Ekebergia_capensis
Elaeagnus_rhamnoides
Emmotum_fagifolium
Enterolobium_schomburgkii
Epilobium_palustre
Equisetum_arvense
Equisetum_fluviatile
Equisetum_palustre
Erica_tetralix
Erythrophleum_chlorostachys
Eschweilera_decolorans
Eschweilera_parviflora
Eschweilera_sagotiana
Eucalyptus_acmenoides
Eucalyptus_grandis
Eucalyptus_umbra
Euphorbia_brachycera
Euphorbia_corollata
Faramea_occidentalis
Fragaria_vesca
Frangula_californica
Fraxinus_uhdei
Galium_verum
Garcinia_macrophylla
Grewia_monticola
Griselinia_littoralis
Guatteria_dumetorum
Gymnosporia_heterophylla
Hakea_cygna
Hakea_preissii
Hamamelis_virginiana
Hampea_appendiculata
Hampea_nutricia
Heliocarpus_appendiculatus
Heracleum_sphondylium
Heteromeles_arbutifolia
Holcus_mollis
Holodiscus_discolor
Ilex_mitis
Inga_acreana
Inga_goldmanii
Inga_nobilis
Inga_sapindoides
Iris_missouriensis
Iryanthera_juruensis
Jacaranda_copaia
Jacaratia_digitata
Jacobaea_vulgaris
Juncus_articulatus
Juniperus_communis
Kigelia_africana
Koeleria_macrantha
Lacmellea_panamensis
Lactuca_muralis
Laetia_corymbulosa
Larix_kaempferi
Larrea_cuneifolia
Larrea_divaricata
Lathyrus_pratensis
Lecythis_zabucajo
Leiospermum_racemosum
Lepechinia_calycina
Leptospermum_polygalifolium
Leucadendron_salignum
Ligustrum_lucidum
Liquidambar_styraciflua
Lotus_pedunculatus
Lysimachia_vulgaris
Lythrum_salicaria
Macaranga_hypoleuca
Machilus_chinensis
Macrolobium_gracile
Magnolia_grandiflora
Malus_sylvestris
Mangifera_indica
Margaritaria_nobilis
Melaleuca_leucadendra
Metrosideros_polymorpha
Mimulus_aurantiacus
Molinia_caerulea
Mosannona_garwoodii
Murraya_paniculata
Myrcia_splendens
Nothofagus_betuloides
Nothofagus_dombeyi
Nothofagus_fusca
Nothofagus_solandri
Ochroma_pyramidale
Ocotea_floribunda
Ocotea_puberula
Opuntia_sulphurea
Ormosia_macrocalyx
Orthion_oblanceolatum
Ouratea_lucens
Oxalis_acetosella
Pachira_aquatica
Pachira_insignis
Packera_multilobata
Palicourea_tetragona
Parinari_campestris
Parkinsonia_praecox
Parnassia_palustris
Peltogyne_venosa
Perebea_xanthochyma
Petrea_volubilis
Picea_sitchensis
Picramnia_latifolia
Pimpinella_major
Pinus_clausa
Pinus_koraiensis
Pinus_massoniana
Piper_aequale
Piscidia_carthagenensis
Pistacia_lentiscus
Pongamia_pinnata
Populus_nigra
Populus_tremuloides
Pourouma_bicolor
Pouteria_cladantha
Pouteria_ephedrantha
Pradosia_cochlearia
Prosopis_flexuosa
Prosopis_torquata
Protium_tenuifolium
Prunus_mahaleb
Prunus_pensylvanica
Pseudocymopterus_montanus
Psidium_guajava
Psychotria_marginata
Pterocarpus_rohrii
Pterygota_alata
Quararibea_asterolepis
Quassia_amara
Quercus_lusitanica
Quercus_robur
Quercus_shumardii
Quercus_stellata
Quercus_variabilis
Ranunculus_bulbosus
Reynoutria_japonica
Rhamnus_alaternus
Rollinia_pittieri
Rosa_woodsii
Rubus_parviflorus
Ruellia_humilis
Rumex_obtusifolius
Salix_cinerea
Salix_lasiolepis
Salix_nigra
Schinopsis_marginata
Schizolobium_parahyba
Sequoia_sempervirens
Shorea_leprosula
Shorea_macroptera
Shorea_parvifolia
Shorea_smithiana
Silphium_terebinthinaceum
Sloanea_guianensis
Sloanea_terniflora
Solanum_dulcamara
Sorbus_aria
Sorocea_steinbachii
Spondias_mombin
Spondias_pinnata
Sporobolus_interruptus
Sterculia_speciosa
Stipa_comata
Tabernaemontana_arborea
Tabernaemontana_donnell-smithii
Tachigali_versicolor
Talisia_princeps
Tanacetum_corymbosum
Tapirira_obtusa
Taraxacum_campylodes
Taxodium_distichum
Taxus_baccata
Terminalia_amazonia
Theobroma_cacao
Theobroma_speciosum
Thevetia_ahouai
Trattinnickia_aspera
Triadica_cochinchinensis
Trichilia_pleeana
Trichilia_quadrijuga
Trichospermum_mexicanum
Trifolium_dubium
Triplaris_cumingiana
Ulmus_glabra
Ulmus_minor
Ulmus_pumila
Ulmus_rubra
Unonopsis_rufescens
Vaccinium_myrtillus
Vaccinium_vitis-idaea
Verbena_bracteata
Veronica_chamaedrys
Viburnum_opulus
Victoria_amazonica
Viola_canina
Vismia_baccifera
Vitex_negundo
Vitex_pinnata
Ziziphus_mucronata
Zuelania_guidonia
GEIGER-fitted comparative model of continuous data
fitted ‘BM’ model parameters:
sigsq = 0.007245
z0 = -1.725638
model summary:
log-likelihood = -147.272308
AIC = 298.544617
AICc = 298.599411
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(EB.LDMC <- fitContinuous(phylo_traits$phy, LDMC[!is.na(LDMC)], model = "EB") )
Warning: The following tips were not found in 'data' and were dropped from 'phy':
Abies_grandis
Acacia_aroma
Acacia_auriculiformis
Acacia_doratoxylon
Acacia_karroo
Acacia_mearnsii
Acacia_praecox
Acer_monspessulanum
Acer_saccharinum
Acer_saccharum
Albizia_saman
Alliaria_petiolata
Alnus_rhombifolia
Ambrosia_psilostachya
Anagallis_tenella
Andira_inermis
Anemone_nemorosa
Annona_spraguei
Apeiba_tibourbou
Araucaria_araucana
Arenaria_serpyllifolia
Artemisia_californica
Artemisia_dracunculus
Aspidosperma_album
Banksia_marginata
Betula_papyrifera
Betula_platyphylla
Bombax_ceiba
Briza_media
Bromus_hordeaceus
Bromus_sterilis
Brosimum_utile
Brownea_grandiceps
Butea_monosperma
Byrsonima_crassifolia
Calamagrostis_epigejos
Calophyllum_longifolium
Campanula_rotundifolia
Capsella_bursa-pastoris
Carex_binervis
Carex_hirta
Carex_pallescens
Carex_panicea
Carya_glabra
Caryocar_glabrum
Casearia_sylvestris
Castanea_dentata
Ceanothus_spinosus
Cecropia_obtusifolia
Cecropia_peltata
Ceiba_speciosa
Cercocarpus_montanus_var._glaber
Cinnamomum_subavenium
Cirsium_wheeleri
Cojoba_rufescens
Colophospermum_mopane
Cordia_bicolor
Cordia_caffra
Cordia_megalantha
Cornus_glabrata
Cornus_sericea
Crepis_pyrenaica
Croton_capitatus
Cupania_rufescens
Cupressus_sempervirens
Cussonia_spicata
Cytisus_scoparius
Dacrydium_cupressinum
Danthonia_decumbens
Daphniphyllum_pentandrum
Drimys_winteri
Dryandra_sessilis
Drypetes_standleyi
Duguetia_quitarensis
Ekebergia_capensis
Elaeagnus_rhamnoides
Emmotum_fagifolium
Enterolobium_schomburgkii
Epilobium_palustre
Equisetum_arvense
Equisetum_fluviatile
Equisetum_palustre
Erica_tetralix
Erythrophleum_chlorostachys
Eschweilera_decolorans
Eschweilera_parviflora
Eschweilera_sagotiana
Eucalyptus_acmenoides
Eucalyptus_grandis
Eucalyptus_umbra
Euphorbia_brachycera
Euphorbia_corollata
Faramea_occidentalis
Fragaria_vesca
Frangula_californica
Fraxinus_uhdei
Galium_verum
Garcinia_macrophylla
Grewia_monticola
Griselinia_littoralis
Guatteria_dumetorum
Gymnosporia_heterophylla
Hakea_cygna
Hakea_preissii
Hamamelis_virginiana
Hampea_appendiculata
Hampea_nutricia
Heliocarpus_appendiculatus
Heracleum_sphondylium
Heteromeles_arbutifolia
Holcus_mollis
Holodiscus_discolor
Ilex_mitis
Inga_acreana
Inga_goldmanii
Inga_nobilis
Inga_sapindoides
Iris_missouriensis
Iryanthera_juruensis
Jacaranda_copaia
Jacaratia_digitata
Jacobaea_vulgaris
Juncus_articulatus
Juniperus_communis
Kigelia_africana
Koeleria_macrantha
Lacmellea_panamensis
Lactuca_muralis
Laetia_corymbulosa
Larix_kaempferi
Larrea_cuneifolia
Larrea_divaricata
Lathyrus_pratensis
Lecythis_zabucajo
Leiospermum_racemosum
Lepechinia_calycina
Leptospermum_polygalifolium
Leucadendron_salignum
Ligustrum_lucidum
Liquidambar_styraciflua
Lotus_pedunculatus
Lysimachia_vulgaris
Lythrum_salicaria
Macaranga_hypoleuca
Machilus_chinensis
Macrolobium_gracile
Magnolia_grandiflora
Malus_sylvestris
Mangifera_indica
Margaritaria_nobilis
Melaleuca_leucadendra
Metrosideros_polymorpha
Mimulus_aurantiacus
Molinia_caerulea
Mosannona_garwoodii
Murraya_paniculata
Myrcia_splendens
Nothofagus_betuloides
Nothofagus_dombeyi
Nothofagus_fusca
Nothofagus_solandri
Ochroma_pyramidale
Ocotea_floribunda
Ocotea_puberula
Opuntia_sulphurea
Ormosia_macrocalyx
Orthion_oblanceolatum
Ouratea_lucens
Oxalis_acetosella
Pachira_aquatica
Pachira_insignis
Packera_multilobata
Palicourea_tetragona
Parinari_campestris
Parkinsonia_praecox
Parnassia_palustris
Peltogyne_venosa
Perebea_xanthochyma
Petrea_volubilis
Picea_sitchensis
Picramnia_latifolia
Pimpinella_major
Pinus_clausa
Pinus_koraiensis
Pinus_massoniana
Piper_aequale
Piscidia_carthagenensis
Pistacia_lentiscus
Pongamia_pinnata
Populus_nigra
Populus_tremuloides
Pourouma_bicolor
Pouteria_cladantha
Pouteria_ephedrantha
Pradosia_cochlearia
Prosopis_flexuosa
Prosopis_torquata
Protium_tenuifolium
Prunus_mahaleb
Prunus_pensylvanica
Pseudocymopterus_montanus
Psidium_guajava
Psychotria_marginata
Pterocarpus_rohrii
Pterygota_alata
Quararibea_asterolepis
Quassia_amara
Quercus_lusitanica
Quercus_robur
Quercus_shumardii
Quercus_stellata
Quercus_variabilis
Ranunculus_bulbosus
Reynoutria_japonica
Rhamnus_alaternus
Rollinia_pittieri
Rosa_woodsii
Rubus_parviflorus
Ruellia_humilis
Rumex_obtusifolius
Salix_cinerea
Salix_lasiolepis
Salix_nigra
Schinopsis_marginata
Schizolobium_parahyba
Sequoia_sempervirens
Shorea_leprosula
Shorea_macroptera
Shorea_parvifolia
Shorea_smithiana
Silphium_terebinthinaceum
Sloanea_guianensis
Sloanea_terniflora
Solanum_dulcamara
Sorbus_aria
Sorocea_steinbachii
Spondias_mombin
Spondias_pinnata
Sporobolus_interruptus
Sterculia_speciosa
Stipa_comata
Tabernaemontana_arborea
Tabernaemontana_donnell-smithii
Tachigali_versicolor
Talisia_princeps
Tanacetum_corymbosum
Tapirira_obtusa
Taraxacum_campylodes
Taxodium_distichum
Taxus_baccata
Terminalia_amazonia
Theobroma_cacao
Theobroma_speciosum
Thevetia_ahouai
Trattinnickia_aspera
Triadica_cochinchinensis
Trichilia_pleeana
Trichilia_quadrijuga
Trichospermum_mexicanum
Trifolium_dubium
Triplaris_cumingiana
Ulmus_glabra
Ulmus_minor
Ulmus_pumila
Ulmus_rubra
Unonopsis_rufescens
Vaccinium_myrtillus
Vaccinium_vitis-idaea
Verbena_bracteata
Veronica_chamaedrys
Viburnum_opulus
Victoria_amazonica
Viola_canina
Vismia_baccifera
Vitex_negundo
Vitex_pinnata
Ziziphus_mucronata
Zuelania_guidoniaWarning:
Parameter estimates appear at bounds:
a
GEIGER-fitted comparative model of continuous data
fitted ‘EB’ model parameters:
a = -0.000001
sigsq = 0.007246
z0 = -1.725640
model summary:
log-likelihood = -147.274172
AIC = 300.548344
AICc = 300.658436
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 39
frequency of best fit = 0.390
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(OU.LDMC <- fitContinuous(phylo_traits$phy, LDMC[!is.na(LDMC)], model = "OU") )
Warning: The following tips were not found in 'data' and were dropped from 'phy':
Abies_grandis
Acacia_aroma
Acacia_auriculiformis
Acacia_doratoxylon
Acacia_karroo
Acacia_mearnsii
Acacia_praecox
Acer_monspessulanum
Acer_saccharinum
Acer_saccharum
Albizia_saman
Alliaria_petiolata
Alnus_rhombifolia
Ambrosia_psilostachya
Anagallis_tenella
Andira_inermis
Anemone_nemorosa
Annona_spraguei
Apeiba_tibourbou
Araucaria_araucana
Arenaria_serpyllifolia
Artemisia_californica
Artemisia_dracunculus
Aspidosperma_album
Banksia_marginata
Betula_papyrifera
Betula_platyphylla
Bombax_ceiba
Briza_media
Bromus_hordeaceus
Bromus_sterilis
Brosimum_utile
Brownea_grandiceps
Butea_monosperma
Byrsonima_crassifolia
Calamagrostis_epigejos
Calophyllum_longifolium
Campanula_rotundifolia
Capsella_bursa-pastoris
Carex_binervis
Carex_hirta
Carex_pallescens
Carex_panicea
Carya_glabra
Caryocar_glabrum
Casearia_sylvestris
Castanea_dentata
Ceanothus_spinosus
Cecropia_obtusifolia
Cecropia_peltata
Ceiba_speciosa
Cercocarpus_montanus_var._glaber
Cinnamomum_subavenium
Cirsium_wheeleri
Cojoba_rufescens
Colophospermum_mopane
Cordia_bicolor
Cordia_caffra
Cordia_megalantha
Cornus_glabrata
Cornus_sericea
Crepis_pyrenaica
Croton_capitatus
Cupania_rufescens
Cupressus_sempervirens
Cussonia_spicata
Cytisus_scoparius
Dacrydium_cupressinum
Danthonia_decumbens
Daphniphyllum_pentandrum
Drimys_winteri
Dryandra_sessilis
Drypetes_standleyi
Duguetia_quitarensis
Ekebergia_capensis
Elaeagnus_rhamnoides
Emmotum_fagifolium
Enterolobium_schomburgkii
Epilobium_palustre
Equisetum_arvense
Equisetum_fluviatile
Equisetum_palustre
Erica_tetralix
Erythrophleum_chlorostachys
Eschweilera_decolorans
Eschweilera_parviflora
Eschweilera_sagotiana
Eucalyptus_acmenoides
Eucalyptus_grandis
Eucalyptus_umbra
Euphorbia_brachycera
Euphorbia_corollata
Faramea_occidentalis
Fragaria_vesca
Frangula_californica
Fraxinus_uhdei
Galium_verum
Garcinia_macrophylla
Grewia_monticola
Griselinia_littoralis
Guatteria_dumetorum
Gymnosporia_heterophylla
Hakea_cygna
Hakea_preissii
Hamamelis_virginiana
Hampea_appendiculata
Hampea_nutricia
Heliocarpus_appendiculatus
Heracleum_sphondylium
Heteromeles_arbutifolia
Holcus_mollis
Holodiscus_discolor
Ilex_mitis
Inga_acreana
Inga_goldmanii
Inga_nobilis
Inga_sapindoides
Iris_missouriensis
Iryanthera_juruensis
Jacaranda_copaia
Jacaratia_digitata
Jacobaea_vulgaris
Juncus_articulatus
Juniperus_communis
Kigelia_africana
Koeleria_macrantha
Lacmellea_panamensis
Lactuca_muralis
Laetia_corymbulosa
Larix_kaempferi
Larrea_cuneifolia
Larrea_divaricata
Lathyrus_pratensis
Lecythis_zabucajo
Leiospermum_racemosum
Lepechinia_calycina
Leptospermum_polygalifolium
Leucadendron_salignum
Ligustrum_lucidum
Liquidambar_styraciflua
Lotus_pedunculatus
Lysimachia_vulgaris
Lythrum_salicaria
Macaranga_hypoleuca
Machilus_chinensis
Macrolobium_gracile
Magnolia_grandiflora
Malus_sylvestris
Mangifera_indica
Margaritaria_nobilis
Melaleuca_leucadendra
Metrosideros_polymorpha
Mimulus_aurantiacus
Molinia_caerulea
Mosannona_garwoodii
Murraya_paniculata
Myrcia_splendens
Nothofagus_betuloides
Nothofagus_dombeyi
Nothofagus_fusca
Nothofagus_solandri
Ochroma_pyramidale
Ocotea_floribunda
Ocotea_puberula
Opuntia_sulphurea
Ormosia_macrocalyx
Orthion_oblanceolatum
Ouratea_lucens
Oxalis_acetosella
Pachira_aquatica
Pachira_insignis
Packera_multilobata
Palicourea_tetragona
Parinari_campestris
Parkinsonia_praecox
Parnassia_palustris
Peltogyne_venosa
Perebea_xanthochyma
Petrea_volubilis
Picea_sitchensis
Picramnia_latifolia
Pimpinella_major
Pinus_clausa
Pinus_koraiensis
Pinus_massoniana
Piper_aequale
Piscidia_carthagenensis
Pistacia_lentiscus
Pongamia_pinnata
Populus_nigra
Populus_tremuloides
Pourouma_bicolor
Pouteria_cladantha
Pouteria_ephedrantha
Pradosia_cochlearia
Prosopis_flexuosa
Prosopis_torquata
Protium_tenuifolium
Prunus_mahaleb
Prunus_pensylvanica
Pseudocymopterus_montanus
Psidium_guajava
Psychotria_marginata
Pterocarpus_rohrii
Pterygota_alata
Quararibea_asterolepis
Quassia_amara
Quercus_lusitanica
Quercus_robur
Quercus_shumardii
Quercus_stellata
Quercus_variabilis
Ranunculus_bulbosus
Reynoutria_japonica
Rhamnus_alaternus
Rollinia_pittieri
Rosa_woodsii
Rubus_parviflorus
Ruellia_humilis
Rumex_obtusifolius
Salix_cinerea
Salix_lasiolepis
Salix_nigra
Schinopsis_marginata
Schizolobium_parahyba
Sequoia_sempervirens
Shorea_leprosula
Shorea_macroptera
Shorea_parvifolia
Shorea_smithiana
Silphium_terebinthinaceum
Sloanea_guianensis
Sloanea_terniflora
Solanum_dulcamara
Sorbus_aria
Sorocea_steinbachii
Spondias_mombin
Spondias_pinnata
Sporobolus_interruptus
Sterculia_speciosa
Stipa_comata
Tabernaemontana_arborea
Tabernaemontana_donnell-smithii
Tachigali_versicolor
Talisia_princeps
Tanacetum_corymbosum
Tapirira_obtusa
Taraxacum_campylodes
Taxodium_distichum
Taxus_baccata
Terminalia_amazonia
Theobroma_cacao
Theobroma_speciosum
Thevetia_ahouai
Trattinnickia_aspera
Triadica_cochinchinensis
Trichilia_pleeana
Trichilia_quadrijuga
Trichospermum_mexicanum
Trifolium_dubium
Triplaris_cumingiana
Ulmus_glabra
Ulmus_minor
Ulmus_pumila
Ulmus_rubra
Unonopsis_rufescens
Vaccinium_myrtillus
Vaccinium_vitis-idaea
Verbena_bracteata
Veronica_chamaedrys
Viburnum_opulus
Victoria_amazonica
Viola_canina
Vismia_baccifera
Vitex_negundo
Vitex_pinnata
Ziziphus_mucronata
Zuelania_guidonia
GEIGER-fitted comparative model of continuous data
fitted ‘OU’ model parameters:
alpha = 0.043364
sigsq = 0.015990
z0 = -1.648994
model summary:
log-likelihood = -94.773735
AIC = 195.547471
AICc = 195.657563
free parameters = 3
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 34
frequency of best fit = 0.340
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
(null.LDMC <- fitContinuous(phylo_traits$phy, LDMC[!is.na(LDMC)], model = "white"))
Warning: The following tips were not found in 'data' and were dropped from 'phy':
Abies_grandis
Acacia_aroma
Acacia_auriculiformis
Acacia_doratoxylon
Acacia_karroo
Acacia_mearnsii
Acacia_praecox
Acer_monspessulanum
Acer_saccharinum
Acer_saccharum
Albizia_saman
Alliaria_petiolata
Alnus_rhombifolia
Ambrosia_psilostachya
Anagallis_tenella
Andira_inermis
Anemone_nemorosa
Annona_spraguei
Apeiba_tibourbou
Araucaria_araucana
Arenaria_serpyllifolia
Artemisia_californica
Artemisia_dracunculus
Aspidosperma_album
Banksia_marginata
Betula_papyrifera
Betula_platyphylla
Bombax_ceiba
Briza_media
Bromus_hordeaceus
Bromus_sterilis
Brosimum_utile
Brownea_grandiceps
Butea_monosperma
Byrsonima_crassifolia
Calamagrostis_epigejos
Calophyllum_longifolium
Campanula_rotundifolia
Capsella_bursa-pastoris
Carex_binervis
Carex_hirta
Carex_pallescens
Carex_panicea
Carya_glabra
Caryocar_glabrum
Casearia_sylvestris
Castanea_dentata
Ceanothus_spinosus
Cecropia_obtusifolia
Cecropia_peltata
Ceiba_speciosa
Cercocarpus_montanus_var._glaber
Cinnamomum_subavenium
Cirsium_wheeleri
Cojoba_rufescens
Colophospermum_mopane
Cordia_bicolor
Cordia_caffra
Cordia_megalantha
Cornus_glabrata
Cornus_sericea
Crepis_pyrenaica
Croton_capitatus
Cupania_rufescens
Cupressus_sempervirens
Cussonia_spicata
Cytisus_scoparius
Dacrydium_cupressinum
Danthonia_decumbens
Daphniphyllum_pentandrum
Drimys_winteri
Dryandra_sessilis
Drypetes_standleyi
Duguetia_quitarensis
Ekebergia_capensis
Elaeagnus_rhamnoides
Emmotum_fagifolium
Enterolobium_schomburgkii
Epilobium_palustre
Equisetum_arvense
Equisetum_fluviatile
Equisetum_palustre
Erica_tetralix
Erythrophleum_chlorostachys
Eschweilera_decolorans
Eschweilera_parviflora
Eschweilera_sagotiana
Eucalyptus_acmenoides
Eucalyptus_grandis
Eucalyptus_umbra
Euphorbia_brachycera
Euphorbia_corollata
Faramea_occidentalis
Fragaria_vesca
Frangula_californica
Fraxinus_uhdei
Galium_verum
Garcinia_macrophylla
Grewia_monticola
Griselinia_littoralis
Guatteria_dumetorum
Gymnosporia_heterophylla
Hakea_cygna
Hakea_preissii
Hamamelis_virginiana
Hampea_appendiculata
Hampea_nutricia
Heliocarpus_appendiculatus
Heracleum_sphondylium
Heteromeles_arbutifolia
Holcus_mollis
Holodiscus_discolor
Ilex_mitis
Inga_acreana
Inga_goldmanii
Inga_nobilis
Inga_sapindoides
Iris_missouriensis
Iryanthera_juruensis
Jacaranda_copaia
Jacaratia_digitata
Jacobaea_vulgaris
Juncus_articulatus
Juniperus_communis
Kigelia_africana
Koeleria_macrantha
Lacmellea_panamensis
Lactuca_muralis
Laetia_corymbulosa
Larix_kaempferi
Larrea_cuneifolia
Larrea_divaricata
Lathyrus_pratensis
Lecythis_zabucajo
Leiospermum_racemosum
Lepechinia_calycina
Leptospermum_polygalifolium
Leucadendron_salignum
Ligustrum_lucidum
Liquidambar_styraciflua
Lotus_pedunculatus
Lysimachia_vulgaris
Lythrum_salicaria
Macaranga_hypoleuca
Machilus_chinensis
Macrolobium_gracile
Magnolia_grandiflora
Malus_sylvestris
Mangifera_indica
Margaritaria_nobilis
Melaleuca_leucadendra
Metrosideros_polymorpha
Mimulus_aurantiacus
Molinia_caerulea
Mosannona_garwoodii
Murraya_paniculata
Myrcia_splendens
Nothofagus_betuloides
Nothofagus_dombeyi
Nothofagus_fusca
Nothofagus_solandri
Ochroma_pyramidale
Ocotea_floribunda
Ocotea_puberula
Opuntia_sulphurea
Ormosia_macrocalyx
Orthion_oblanceolatum
Ouratea_lucens
Oxalis_acetosella
Pachira_aquatica
Pachira_insignis
Packera_multilobata
Palicourea_tetragona
Parinari_campestris
Parkinsonia_praecox
Parnassia_palustris
Peltogyne_venosa
Perebea_xanthochyma
Petrea_volubilis
Picea_sitchensis
Picramnia_latifolia
Pimpinella_major
Pinus_clausa
Pinus_koraiensis
Pinus_massoniana
Piper_aequale
Piscidia_carthagenensis
Pistacia_lentiscus
Pongamia_pinnata
Populus_nigra
Populus_tremuloides
Pourouma_bicolor
Pouteria_cladantha
Pouteria_ephedrantha
Pradosia_cochlearia
Prosopis_flexuosa
Prosopis_torquata
Protium_tenuifolium
Prunus_mahaleb
Prunus_pensylvanica
Pseudocymopterus_montanus
Psidium_guajava
Psychotria_marginata
Pterocarpus_rohrii
Pterygota_alata
Quararibea_asterolepis
Quassia_amara
Quercus_lusitanica
Quercus_robur
Quercus_shumardii
Quercus_stellata
Quercus_variabilis
Ranunculus_bulbosus
Reynoutria_japonica
Rhamnus_alaternus
Rollinia_pittieri
Rosa_woodsii
Rubus_parviflorus
Ruellia_humilis
Rumex_obtusifolius
Salix_cinerea
Salix_lasiolepis
Salix_nigra
Schinopsis_marginata
Schizolobium_parahyba
Sequoia_sempervirens
Shorea_leprosula
Shorea_macroptera
Shorea_parvifolia
Shorea_smithiana
Silphium_terebinthinaceum
Sloanea_guianensis
Sloanea_terniflora
Solanum_dulcamara
Sorbus_aria
Sorocea_steinbachii
Spondias_mombin
Spondias_pinnata
Sporobolus_interruptus
Sterculia_speciosa
Stipa_comata
Tabernaemontana_arborea
Tabernaemontana_donnell-smithii
Tachigali_versicolor
Talisia_princeps
Tanacetum_corymbosum
Tapirira_obtusa
Taraxacum_campylodes
Taxodium_distichum
Taxus_baccata
Terminalia_amazonia
Theobroma_cacao
Theobroma_speciosum
Thevetia_ahouai
Trattinnickia_aspera
Triadica_cochinchinensis
Trichilia_pleeana
Trichilia_quadrijuga
Trichospermum_mexicanum
Trifolium_dubium
Triplaris_cumingiana
Ulmus_glabra
Ulmus_minor
Ulmus_pumila
Ulmus_rubra
Unonopsis_rufescens
Vaccinium_myrtillus
Vaccinium_vitis-idaea
Verbena_bracteata
Veronica_chamaedrys
Viburnum_opulus
Victoria_amazonica
Viola_canina
Vismia_baccifera
Vitex_negundo
Vitex_pinnata
Ziziphus_mucronata
Zuelania_guidonia
GEIGER-fitted comparative model of continuous data
fitted ‘white’ model parameters:
sigsq = 0.182623
z0 = -1.591252
model summary:
log-likelihood = -126.267505
AIC = 256.535011
AICc = 256.589805
free parameters = 2
Convergence diagnostics:
optimization iterations = 100
failed iterations = 0
number of iterations with same best fit = 100
frequency of best fit = 1.000
object summary:
'lik' -- likelihood function
'bnd' -- bounds for likelihood search
'res' -- optimization iteration summary
'opt' -- maximum likelihood parameter estimates
Moderate signal
Check AIC weights
(aics_models <- setNames(c(AIC(BM.Height), AIC(EB.Height), AIC(OU.Height), AIC(null.Height)), c("BM", "EB", "OU", "NULL")) )
BM EB OU NULL
1901.959 1903.966 1823.619 2160.758
aic.w(aics_models)
BM EB OU NULL
0 0 1 0
Again, OU
So overall , the Ornstein-Uhlenbeck model seems to be the best model for most of our data, except Nmass which varies along a grand mean (white-noise).
Check the categorical trait Growth Form
We will use the \(\delta\) statistic for this.
source("https://raw.githubusercontent.com/mrborges23/delta_statistic/master/code.R")
trait <- filtered_data %>% dplyr::select(Species, `Growth Form`) %>% column_to_rownames("Species") %>% .[rownames(phylo_traits$data),]
# This takes a very long time, go for lunch or something ...
#deltaA <- delta(trait,phylo_traits$phy,0.1,0.0589,10000,10,100)
And check the p-value by using randomization.
This can take a while again …
# random_delta <- rep(NA,100)
# for (i in 1:100){
# rtrait <- sample(trait)
# random_delta[i] <- delta(rtrait,phylo_traits$phy,0.1,0.0589,10000,10,100)
# }
# p_value <- sum(random_delta>deltaA)/length(random_delta)
# boxplot(random_delta)
# abline(h=deltaA,col="red")
For binary data you can uncomment the code bellow
So now that we have an idea of the strength of the phylogenetic signal, it’s time to decide on a method of imputation.
Currently there are 3 main flavours of phylogenetic imputation see review in Molna-Venegas et al. 2018:
pGLM can easily be implemented by fitting a pVCV matrix as the
correlation structure in a gls model (package nlme). This
assumes a Brownian model of evolution. More sophisticated ways of
implementing a pGLM can be done with the package rPhylopars
and several evolutionary models can be fitted. See tutorial here
PVR can be fit using a randomForest approach in the package
missForest. One simply needs to transform the phylogeny
into a phylogenetic distance matrix and use this in a PCoA. Then the
PCoA axes are added as extra variables to the trait matrix as
predictors.
PEM can be fit using the R package MPSEM. This expands
on the PVR method by adding an evolutionary model and is thus considered
more appropriate for evo-eco modelling. The tree is transformed prior to
transforming into distance matrix and PCoA.
A fourth option would be to use Bayesian Hierarchical Matrix Factorization (BHPMF), which has been used by many publications for trait imputation. This method however doesn’t really use the phylogeny or quantify phylogenetic signal. It uses a taxonomic hierarchy and assumes that all ranks are equal and contain information. In reality, all taxonomic ranks are not equal, as some genera can be older than certain families for example. As with all Bayesian models, it will usually be much slower but will provide uncertainty of the prediction.
Let’s see a simple example with pGLM
First, let’s randomly remove 10% of all values in the data frame. The
missForest package can make this easy for us using the
prodNA function.
set.seed(4312)
missing <- phylo_traits$data[,-1] %>% prodNA(0.1)
head(missing)
RPhylopars is a bit fussy and requires a
species column (all lower case) that is the very first
column in the trait matrix.
missing <- missing %>%
mutate(species = rownames(.)) %>%
dplyr::select(species, everything())
model.BM <- phylopars(tree = phylo_traits$phy, trait_data = missing, model = "BM")
warning: solve(): system is singular (rcond: 5.60546e-27); attempting approx solution
warning: solve(): system is singular (rcond: 1.80444e-33); attempting approx solution
warning: solve(): system is singular (rcond: 1.54393e-32); attempting approx solution
warning: solve(): system is singular (rcond: 3.53749e-29); attempting approx solution
model.BM
Phylogenetic trait variance-covariance
LA Nmass LMA SSD Height Seedmass
LA 0.146855292 -0.002810604 -0.005329487 -0.004416941 0.044477787 0.052505321
Nmass -0.002810604 0.060812588 -0.019895646 -0.018075167 0.033858310 0.029082235
LMA -0.005329487 -0.019895646 0.021178855 0.008857747 -0.013069464 -0.001161069
SSD -0.004416941 -0.018075167 0.008857747 0.012493858 -0.003575309 -0.005736158
Height 0.044477787 0.033858310 -0.013069464 -0.003575309 0.080096359 0.042377476
Seedmass 0.052505321 0.029082235 -0.001161069 -0.005736158 0.042377476 0.139364874
Brownian motion model
summary(model.BM)
Brownian motion model
phylogenetic mean phylogenetic sd
LA 5.3937552 0.3832
Nmass 2.9001291 0.2466
LMA 4.5167510 0.1455
SSD -1.2674510 0.1118
Height 0.7516388 0.2830
Seedmass -2.6710471 0.3733
We can compare with other models. EB, OU, lambda and two variants of multivariate-OU (with full alpha or fixed alpha which corresponds to the adaptation rate parameter)
model.OU <- phylopars(trait_data = missing, tree = phylo_traits$phy, model = "OU")
model.mvOU <- phylopars(trait_data = missing, tree = phylo_traits$phy, model = "mvOU", full_alpha = TRUE, usezscores = FALSE) # the usezscores = F is a short term solution to a version issues on CRAN, see here https://github.com/ericgoolsby/Rphylopars/issues/54
model.mvOU_diag <- phylopars(trait_data = missing, tree = phylo_traits$phy, model = "mvOU", full_alpha = FALSE, usezscores = FALSE)
warning: solve(): system is singular (rcond: 2.35824e-132); attempting approx solution
warning: solve(): system is singular (rcond: 6.37787e-189); attempting approx solution
warning: solve(): system is singular (rcond: 2.25478e-23); attempting approx solution
warning: solve(): system is singular (rcond: 9.93784e-24); attempting approx solution
warning: solve(): system is singular (rcond: 1.88375e-206); attempting approx solution
warning: solve(): system is singular (rcond: 1.88375e-206); attempting approx solution
warning: solve(): system is singular (rcond: 5.16434e-234); attempting approx solution
warning: solve(): system is singular (rcond: 5.16434e-234); attempting approx solution
warning: solve(): system is singular (rcond: 6.6632e-32); attempting approx solution
warning: solve(): system is singular (rcond: 6.6632e-32); attempting approx solution
warning: solve(): system is singular (rcond: 1.83029e-18); attempting approx solution
warning: solve(): system is singular (rcond: 1.83029e-18); attempting approx solution
warning: solve(): system is singular (rcond: 2.59378e-57); attempting approx solution
warning: solve(): system is singular (rcond: 2.59378e-57); attempting approx solution
warning: solve(): system is singular (rcond: 2.59431e-57); attempting approx solution
warning: solve(): system is singular (rcond: 5.84998e-71); attempting approx solution
warning: solve(): system is singular (rcond: 5.84998e-71); attempting approx solution
warning: solve(): system is singular (rcond: 5.84998e-71); attempting approx solution
warning: solve(): system is singular (rcond: 4.97832e-18); attempting approx solution
warning: solve(): system is singular (rcond: 4.97832e-18); attempting approx solution
warning: solve(): system is singular (rcond: 4.97832e-18); attempting approx solution
warning: solve(): system is singular (rcond: 8.83805e-26); attempting approx solution
warning: solve(): system is singular (rcond: 8.83805e-26); attempting approx solution
warning: solve(): system is singular (rcond: 1.56371e-21); attempting approx solution
warning: solve(): system is singular (rcond: 1.56537e-21); attempting approx solution
model.EB <- phylopars(trait_data = missing, tree = phylo_traits$phy, model = "EB")
warning: solve(): system is singular (rcond: 4.3832e-17); attempting approx solution
warning: solve(): system is singular (rcond: nan); attempting approx solution
warning: solve(): system is singular (rcond: 1.38784e-39); attempting approx solution
warning: solve(): system is singular (rcond: nan); attempting approx solution
warning: solve(): system is singular (rcond: 2.30631e-27); attempting approx solution
warning: solve(): system is singular (rcond: nan); attempting approx solution
warning: solve(): system is singular (rcond: 5.0265e-108); attempting approx solution
warning: solve(): system is singular (rcond: nan); attempting approx solution
warning: solve(): system is singular (rcond: 3.43122e-54); attempting approx solution
warning: solve(): system is singular (rcond: 5.41758e-157); attempting approx solution
warning: solve(): system is singular (rcond: 2.66719e-17); attempting approx solution
warning: solve(): system is singular (rcond: 6.56869e-18); attempting approx solution
warning: solve(): system is singular (rcond: 2.42022e-37); attempting approx solution
warning: solve(): system is singular (rcond: 5.4724e-25); attempting approx solution
model.lambda <- phylopars(trait_data = missing, tree = phylo_traits$phy, model = "lambda")
warning: solve(): system is singular (rcond: 1.46268e-85); attempting approx solution
warning: solve(): system is singular (rcond: 3.32532e-215); attempting approx solution
warning: solve(): system is singular (rcond: 1.84798e-25); attempting approx solution
warning: solve(): system is singular (rcond: 1.80378e-65); attempting approx solution
warning: solve(): system is singular (rcond: 3.53515e-18); attempting approx solution
warning: solve(): system is singular (rcond: 7.05819e-87); attempting approx solution
warning: solve(): system is singular (rcond: 9.35436e-22); attempting approx solution
warning: solve(): system is singular (rcond: 4.40538e-72); attempting approx solution
warning: solve(): system is singular (rcond: 1.04763e-31); attempting approx solution
warning: solve(): system is singular (rcond: 6.86129e-42); attempting approx solution
warning: solve(): system is singular (rcond: 1.70055e-19); attempting approx solution
warning: solve(): system is singular (rcond: 1.53256e-18); attempting approx solution
warning: solve(): system is singular (rcond: 6.57095e-28); attempting approx solution
warning: solve(): system is singular (rcond: 6.59432e-27); attempting approx solution
warning: solve(): system is singular (rcond: 6.07831e-78); attempting approx solution
warning: solve(): system is singular (rcond: 6.01054e-17); attempting approx solution
warning: solve(): system is singular (rcond: 4.50702e-33); attempting approx solution
This can take quite a while to run.
We can compare the AICs of each model.
aics_models2 <- setNames(c(AIC(model.BM), AIC(model.OU), AIC(model.mvOU), AIC(model.mvOU_diag), AIC(model.EB), AIC(model.lambda)), c("BM", "OU", "mvOU", "mvOU_diag", "EB", "lambda"))
aics_models2
BM OU mvOU mvOU_diag EB lambda
8502.616 1029265.329 10272.851 6794.119 8505.241 6174.510
aicw(aics_models2)
NA
Here a lambda transformation applied to the tree seems to fit the data best.
We can now get the imputed values from the model and use the variance to figure out their lower and upper 95% conf interval.
model.lambda$anc_recon[1:5,] %>% as_tibble() # means
model.lambda$anc_var[1:5,] %>% as_tibble() # variances
model.lambda$anc_recon[1:5,]- sqrt(model.lambda$anc_var[1:5,])*1.96 # Lower 95% CI
LA Nmass LMA SSD Height Seedmass
Silphium_terebinthinaceum 9.988363 2.667228 5.430707 -1.073132 -2.1754462 2.8608171
Heliopsis_helianthoides 7.715864 3.432373 3.830609 -1.629847 0.2745473 1.5187908
Eupatorium_perfoliatum 9.598143 2.563614 4.274364 -1.483378 -1.6735469 -2.3848803
Gaillardia_pinnatifida 4.794798 2.638274 3.312451 -1.719268 -2.9397255 0.2126082
Hymenoxys_richardsonii 4.941642 2.888147 4.565949 -1.632553 -1.5606477 -4.0153066
model.lambda$anc_recon[1:5,]+ sqrt(model.lambda$anc_var[1:5,])*1.96 #Upper 95% CI
LA Nmass LMA SSD Height Seedmass
Silphium_terebinthinaceum 9.988363 2.667228 5.430707 -1.073132 -2.1754462 2.8608171
Heliopsis_helianthoides 7.715864 3.432373 3.830609 -1.629847 0.2745473 1.5187908
Eupatorium_perfoliatum 9.598143 3.482349 4.274364 -1.483378 1.3419938 -2.3848803
Gaillardia_pinnatifida 4.794798 3.634277 4.746396 -1.719268 -2.9397255 0.2126082
Hymenoxys_richardsonii 4.941642 2.888147 4.565949 -1.632553 -1.5606477 1.3747380
You can see that the variance for Seedmass is quite high. For Solidago_rigida, the trait values go from -2.57 to 3.27.
# variance in Seedmass missing values imputed by Rphylopars
model.lambda$anc_var[which(is.na(missing$Seedmass)),"Seedmass"] %>% summary()
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.589 1.979 2.628 2.881 3.593 5.557
We could use this to filter out imputed values that have too high a variance.
Now we will use the Phylogenetic Eigenvector Regression technique (PVR).
phylo.cor <- cophenetic.phylo(phylo_traits$phy)
# Check if the distance matrix is euclidean, if not, can square root the phylogeny
is.euclid(as.dist(phylo.cor))
[1] TRUE
phylo.pcoa <- dudi.pco(d = as.dist(phylo.cor), scannf = F, full = T)
# from ape, calculates broken stick automatically
phylo.pcoa2 <- pcoa(D = as.dist(phylo.cor))
Ok now we have our phylogenetic eigenevectors, the question is how many do we select to include into our regression?
Some people say all, some people say use a broken stick model, some people say use those that explain at least 50% of variance. See discussion by Rolhf , 2001 and Martins et al. 2002.
For now I will use the broken stick model. Some have suggested
fa.parallel from package psych but I have had
many issues with this. We can check this out visually and extract these
from the object
ggplot() + geom_point(data=phylo.pcoa2$values, aes(x=1:nrow(phylo.pcoa2$values), y=Relative_eig), col="blue") + geom_point(data=phylo.pcoa2$values, aes(x=1:nrow(phylo.pcoa2$values), y=Broken_stick), col="red")
(sig.EV <-phylo.pcoa2$values[which(phylo.pcoa2$values$Relative_eig > phylo.pcoa2$values$Broken_stick),])
# Take the eigenvectors and add them to the traits matrix
missing_phylo <- cbind(missing, phylo.pcoa$li[,1:nrow(sig.EV)]) %>% dplyr::select(-species)
pvr <- missForest(xmis = missing_phylo, ntree = 1000, maxiter = 100, variablewise = T)
Check the output
head(pvr)
$ximp
$OOBerror
MSE MSE MSE MSE MSE MSE MSE MSE MSE
1.79489167 0.08689951 0.18307350 0.08737272 0.75754404 2.97288929 0.00000000 0.00000000 0.00000000
MSE MSE MSE MSE MSE MSE MSE MSE MSE
0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
MSE MSE MSE MSE MSE MSE MSE MSE MSE
0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
MSE MSE MSE MSE MSE MSE MSE MSE MSE
0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
We can see that the Out Of Bag Error (OOB) is much higher for trait 6 (ie. Seed Mass) with an OOB of 3.03 compared to others.
Compare the results of the true matrix, the rPhylopars approach and the missForest approach
phylo_traits$data[,-1]
pvr$ximp
data.frame(model.lambda$anc_recon[1:200,])
We can check how well this has done for each variable
# LMA
missing.LMA <- which(is.na(missing$LMA))
cor(phylo_traits$data[missing.LMA,"LMA"], pvr$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.84 PVR"
cor(phylo_traits$data[missing.LMA,"LMA"], model.lambda$anc_recon[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
#LA
missing.LA <- which(is.na(missing$LA))
cor(phylo_traits$data[missing.LA,"LA"], pvr$ximp[missing.LA,"LA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.7 PVR"
cor(phylo_traits$data[missing.LA,"LA"], model.lambda$anc_recon[missing.LA,"LA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
Not great here for both of them. pGLM is much better though.
#Nmass
missing.Nmass <- which(is.na(missing$Nmass))
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.6 PVR"
cor(phylo_traits$data[missing.Nmass,"Nmass"], model.lambda$anc_recon[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.64 pGLM"
#SSD
missing.SSD <- which(is.na(missing$SSD))
cor(phylo_traits$data[missing.SSD,"SSD"], pvr$ximp[missing.SSD,"SSD"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.87 PVR"
cor(phylo_traits$data[missing.SSD,"SSD"], model.lambda$anc_recon[missing.SSD,"SSD"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.86 pGLM"
#Seedmass
missing.Seedmass <- which(is.na(missing$Seedmass))
cor(phylo_traits$data[missing.Seedmass,"Seedmass"], pvr$ximp[missing.Seedmass,"Seedmass"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.81 PVR"
cor(phylo_traits$data[missing.Seedmass,"Seedmass"], model.lambda$anc_recon[missing.Seedmass,"Seedmass"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.84 pGLM"
And height
#Height
missing.Height <- which(is.na(missing$Height))
cor(phylo_traits$data[missing.Height,"Height"], pvr$ximp[missing.Height,"Height"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.88 PVR"
cor(phylo_traits$data[missing.Height,"Height"], model.lambda$anc_recon[missing.Height,"Height"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.86 pGLM"
Ok so overall it’s not too bad! Both methods are very comparable and only Leaf Area and Nitrogen Mass are poorly reconstructed.
What happens if we try and impute these without phylogenetic information?
impute.nophylo <- missForest(xmis = missing[,-1], ntree = 1000, maxiter = 100, variablewise = T)
cor(phylo_traits$data[missing.LA,"LA"], impute.nophylo$ximp[missing.LA,"LA"])
[1] 0.4948403
cor(phylo_traits$data[missing.Nmass,"Nmass"], impute.nophylo$ximp[missing.Nmass,"Nmass"])
[1] 0.6032527
Ok so for Leaf Area, it’s a tad worse than when using
phylogeny, for Nmass it’s the same. This is because
Leaf Area doesn’t have super strong phylogenetic signal and
Nmass has no phylogenetic signal at all and is best fitted
by a white-noise model.
It is also possible that the first phylogenetic eigenvectors that we
selected aren’t the ones that best suit Leaf Area, so the
PVR method does poorly.
We can check the out of bag error
rbind(pvr$OOBerror[1:6], impute.nophylo$OOBerror) %>% `colnames<-`(names(impute.nophylo$ximp)) %>% `rownames<-`(c("RF_PVR", "RF_no_phylo"))
LA Nmass LMA SSD Height Seedmass
RF_PVR 1.794892 0.08689951 0.1830735 0.08737272 0.7575440 2.972889
RF_no_phylo 2.106433 0.10523274 0.1649866 0.10091698 0.9563365 5.098242
Interesting, the OOB error for NMass is quite low compared to others, but the correlation of the imputed values is bad.
Let’s quickly check the correlation between traits
cor(phylo_traits$data[,-1])
LA Nmass LMA SSD Height Seedmass
LA 1.00000000 0.09317154 -0.06162798 0.08995956 0.4248158 0.40655217
Nmass 0.09317154 1.00000000 -0.56983896 -0.33316314 -0.1838869 -0.06222018
LMA -0.06162798 -0.56983896 1.00000000 0.56796971 0.3988510 0.31703906
SSD 0.08995956 -0.33316314 0.56796971 1.00000000 0.7194188 0.57530529
Height 0.42481579 -0.18388692 0.39885096 0.71941879 1.0000000 0.67603235
Seedmass 0.40655217 -0.06222018 0.31703906 0.57530529 0.6760323 1.00000000
This part is very much experimental, and I don’t quite have a strong grip on it yet.
library(MPSEM)
Warning: package ‘MPSEM’ was built under R version 4.2.3
tree.pgraph <- Phylo2DirectedGraph(phylo_traits$phy)
tree.PEM <-PEM.build(tree.pgraph, d="distance",sp="species") # here I'm not adding any a or psi parameters
as.data.frame(tree.PEM)
This next solution requires complete trait values, no NAs.
# here we can let the data speak for itself
# use RMLE to estimate a
tree.PEM_opt2 <- PEM.fitSimple(
y = as.matrix(phylo_traits$data[,-1]),
x=NULL,
w = tree.pgraph,
d = "distance", sp="species",
lower = 0, upper = 1)
The package uses a linear model via stepwise selection based on AIC criterion. It can only do it for one trait at a time and can accept a certain amount of NA values.
## Estimating ancestral trait values:
ANCloc <- getAncGraphLocations(tree.pgraph)
PEMfsAnc <- PEM.fitSimple(y=as.matrix(phylo_traits$data[,-1]),
x=NULL,
w=ANCloc$x,
d="distance",sp="species",
lower=0,upper=1)
PEManc1 <- Locations2PEMscores(PEMfsAnc, ANCloc)
# there is another step to this section that's i can't figure how to work out yet
Now use these in the imputation method with RF.
missing_phylo <- cbind(missing, as.data.frame(tree.PEM_opt2)[,1:ncol(sig.EV)]) %>% dplyr::select(-species)
pem <- missForest(xmis = missing_phylo, ntree = 1000, maxiter = 100, variablewise = T)
Check
# LMA
missing.LMA <- which(is.na(missing$LMA))
cor(phylo_traits$data[missing.LMA,"LMA"], pvr$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.84 PVR"
cor(phylo_traits$data[missing.LMA,"LMA"], model.lambda$anc_recon[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
cor(phylo_traits$data[missing.LMA,"LMA"], pem$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.81 PEM"
#LA
missing.LA <- which(is.na(missing$LA))
cor(phylo_traits$data[missing.LA,"LA"], pvr$ximp[missing.LA,"LA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.7 PVR"
cor(phylo_traits$data[missing.LA,"LA"], model.lambda$anc_recon[missing.LA,"LA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
cor(phylo_traits$data[missing.LA,"LA"], pem$ximp[missing.LA,"LA"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.66 PEM"
#Nmass
missing.Nmass <- which(is.na(missing$Nmass))
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.6 PVR"
cor(phylo_traits$data[missing.Nmass,"Nmass"], model.lambda$anc_recon[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.64 pGLM"
cor(phylo_traits$data[missing.Nmass,"Nmass"], pem$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.6 PEM"
So very similar to PVR if not worse, not much of an improvement really … Unless I’m doing this wrong?
One solution is to run several rounds of imputation.
One could use Rphylopars to first impute the continuous variables that have strong phylogenetic signal and then run phylogenetic eigenvectors with RF to impute the remaining categorical traits and traits that have low phylogenetic signal.
Let’s re-run the lambda model but remove NMass. Then we will add Nmass and the phylogenetic eigenvectors.
model.lambda$anc_recon %>% dim()
[1] 999 6
model.lambda2 <- phylopars(trait_data = missing %>% dplyr::select(-Nmass) , tree = phylo_traits$phy, model = "lambda")
warning: solve(): system is singular (rcond: 1.84019e-22); attempting approx solution
warning: solve(): system is singular (rcond: 3.38169e-222); attempting approx solution
warning: solve(): system is singular (rcond: 3.0712e-17); attempting approx solution
warning: solve(): system is singular (rcond: 1.82178e-17); attempting approx solution
warning: solve(): system is singular (rcond: 6.09341e-41); attempting approx solution
warning: solve(): system is singular (rcond: 2.34904e-75); attempting approx solution
warning: solve(): system is singular (rcond: 1.20236e-19); attempting approx solution
warning: solve(): system is singular (rcond: 2.8273e-26); attempting approx solution
warning: solve(): system is singular (rcond: 2.27438e-31); attempting approx solution
Now add the Nmass
comb.model <- data.frame(model.lambda2$anc_recon[1:500,], Nmass=missing$Nmass, phylo.pcoa$li[,1:nrow(sig.EV)])
pvr.comb <- missForest(xmis = comb.model, ntree = 1000, maxiter = 100, variablewise = T)
And now extract these and compare
# LMA
missing.LMA <- which(is.na(missing$LMA))
cor(phylo_traits$data[missing.LMA,"LMA"], pvr$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.84 PVR"
cor(phylo_traits$data[missing.LMA,"LMA"], model.lambda$anc_recon[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
cor(phylo_traits$data[missing.LMA,"LMA"], pem$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.81 PEM"
cor(phylo_traits$data[missing.LMA,"LMA"], pvr.comb$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"combined")
[1] "cor 0.75 combined"
Nmass
#Nmass
missing.Nmass <- which(is.na(missing$Nmass))
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.6 PVR"
cor(phylo_traits$data[missing.Nmass,"Nmass"], model.lambda$anc_recon[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.64 pGLM"
cor(phylo_traits$data[missing.Nmass,"Nmass"], pem$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.6 PEM"
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr.comb$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"combined")
[1] "cor 0.58 combined"
Interesting, it seems like the combined method is actually worse than any method!
What if we include Nmass in the original imputation, but then remove those results and re-run the randomForest ?
model.lambda3 <- phylopars(trait_data = missing , tree = phylo_traits$phy, model = "lambda")
warning: solve(): system is singular (rcond: 1.46268e-85); attempting approx solution
warning: solve(): system is singular (rcond: 3.32532e-215); attempting approx solution
warning: solve(): system is singular (rcond: 1.84798e-25); attempting approx solution
warning: solve(): system is singular (rcond: 1.80378e-65); attempting approx solution
warning: solve(): system is singular (rcond: 3.53515e-18); attempting approx solution
warning: solve(): system is singular (rcond: 7.05819e-87); attempting approx solution
warning: solve(): system is singular (rcond: 9.35436e-22); attempting approx solution
warning: solve(): system is singular (rcond: 4.40538e-72); attempting approx solution
warning: solve(): system is singular (rcond: 1.04763e-31); attempting approx solution
warning: solve(): system is singular (rcond: 6.86129e-42); attempting approx solution
warning: solve(): system is singular (rcond: 1.70055e-19); attempting approx solution
warning: solve(): system is singular (rcond: 1.53256e-18); attempting approx solution
warning: solve(): system is singular (rcond: 6.57095e-28); attempting approx solution
warning: solve(): system is singular (rcond: 6.59432e-27); attempting approx solution
warning: solve(): system is singular (rcond: 6.07831e-78); attempting approx solution
warning: solve(): system is singular (rcond: 6.01054e-17); attempting approx solution
warning: solve(): system is singular (rcond: 4.50702e-33); attempting approx solution
comb.model2 <- data.frame(model.lambda3$anc_recon[1:500,-2], Nmass=missing$Nmass, phylo.pcoa$li[,1:nrow(sig.EV)])
pvr.comb2 <- missForest(xmis = comb.model2, ntree = 1000, maxiter = 100, variablewise = T)
#LMA
cor(phylo_traits$data[missing.LMA,"LMA"], pvr$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.84 PVR"
cor(phylo_traits$data[missing.LMA,"LMA"], model.lambda$anc_recon[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.82 pGLM"
cor(phylo_traits$data[missing.LMA,"LMA"], pem$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.81 PEM"
cor(phylo_traits$data[missing.LMA,"LMA"], pvr.comb2$ximp[missing.LMA,"LMA"]) %>% round(2) %>% paste("cor",.,"combined")
[1] "cor 0.82 combined"
#Nmass
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PVR")
[1] "cor 0.6 PVR"
cor(phylo_traits$data[missing.Nmass,"Nmass"], model.lambda$anc_recon[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"pGLM")
[1] "cor 0.64 pGLM"
cor(phylo_traits$data[missing.Nmass,"Nmass"], pem$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"PEM")
[1] "cor 0.6 PEM"
cor(phylo_traits$data[missing.Nmass,"Nmass"], pvr.comb2$ximp[missing.Nmass,"Nmass"]) %>% round(2) %>% paste("cor",.,"combined")
[1] "cor 0.58 combined"
Ok, so even here it’s not as good. interesting.
Perhaps this is best done only with remaining categorical traits.
Final test
BHPMF.test <- GapFilling(trait.info, hierarchy.info,
mean.gap.filled.output.path = paste0(tmp.dir, "/mean_gap_filled.txt"),
std.gap.filled.output.path = paste0(tmp.dir, "/std_gap_filled.txt"),
tmp.dir = tmp.dir)
./fold1/Tunning
RMSE for the test data: 0.7519486